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Summary 


This is the final report on our one year grant for the study of "Moisture Content .and .Migration 
Dynamics in Unsaturated Porous Media". We proposed fundamental studies of fluid mechanics 
and transport in partially saturated soils. Two studies were envisioned. In the first, an 
undergraduate student would be involved in the solution of transient diffusion problems in support 
of the development of probes for the in-situ measurement of moisture content In the second, 
numerical and analytical methods were to be used to study the fundamental problem of meniscus 
and saturation front propagation in geometric models of porous media. 

The first project involved Mr. James MacDonald, an undergraduate in Chemical Engineering at 
Stanford working closely with Dr. Boris Yendler, an NRC Fellow at NASA Ames. Dr. Yendler 
had responsibility for the day to day progress of Mr. MacDonald: Prof. Homsy reviewed the work 
on a periodic basis and made some technical suggestions regarding the senes expansion and 
inversion of Laplace transforms. Mr. MacDonald and Dr. Yendler have co-authored a paper, a 
copy of which is appended. 

The second project began with the engagement of Mr. Ahmed Farooq, an advanced PhD studentin 
Mechanical Engineering at Stanford. Following some initial discussions between Mr. Farooq, Dr. 
Yendler and Prof. Homsy, Mr. Farooq conducted a short assessment of the applicability of 
standard ’soil science' models to treat the problem of the propagation of satmauon fron^ m porous 
media. Mr. Farooq's progress was not entirely satisfactory, and it was decided that he should 
leave the project His short report is also appended. Following this. Prof. Homsy took the 
opportunity to acquire a flexible finite element program for the solution of meniscuspropagation 1 
Sro^cs. ™? opportunity was afforded by the sabbatical visit of Prof. Bamin Ktomrum from 
Washington University, St. Louis. The program has been acquired and another of Prof. Homsy s 
PhD students, Ms. Jean Ro, has familiarized herself with the program and has completed some 
benchmark runs. A short description of the program and its capabilities is also attached. Althoug 
the budgeted funds have been expended, Ms. Ro is continuing to work on the problemot 
meniscus propagation through geometrically complex pore space, and all pubhcations based on her 
work will carry an appropriate acknowledgement to the NASA grant that allowed the acquisition o 

the finite element program. 
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Project Defin ition and Objectives 

Thiii work was undertaken in response to a problem encountered while determining thermal 
conductivity ch aracteristics of partially saturated soils. Ongoing research utilized a cylindrical “thermal 
probe” containing & heat source and a transducer in aider to measure thermal dissipation within a 
homoge no us pa rtially soil that was packed via protocol (please sec attachment “Moisture 

Measurements" for reasoning of probe selection). An IBM computer was then used to interpret the 
transducer pulses resulting in profiles of temperature versus time. These profiles were then utilized to back- 
the thermal conductivity of the medium. The problem was that the resulting profiles differed 
Henenriiiw on the spatial placement of the sensors. For example, if a sensor wasplaccd m the medium a 
temperature p rofil e would be obtained, but if the same sensor was placed in a different area of die 
h omog ene ous medium, a different profile was obtained. 

The multiple curves show a difference both in magnitude and slope dining short time, but the long 
rim* data only shows a magnitude offset which is presumed to be due to the offset introduced in die short 
time regime. This observed offset is hypothesized to be a result of variations in the contact resistivity m the 
sensor's surface. The variance occurred due to difference in soil packing around the sensor upon each 
insertion This work will attempt to minimize or eliminate die variance in die contact rewshvity effect in 
order to the numerous experimental curves into one accurate curve which will allow the calculation 

of die medium's thermal conductivity. 


Approach To Solution 

Prior to this work, the first attempt to isolate the contact resistivity resulted in applying the 
cylindrical “thermal probe" equations derived by Blackwell 1 to the data. This i nc l uded both the original 
eauadoos and the long time probe and short time medium equations he combined in an atiempt to account 
for the ccnductivityof die probe. Without knowing the heat transfer coefficient, the variable that describes 
die contact resistivity, or the thermal conductivity of the medium, the theoretical equations for short and long 
rim* were solved for both terms. This same procedure was carried out for multiple runs within the same 
medium but failed to reproduce a consistent thermal conductivity. 

This failure brought into question the validity of one of the primary assumptions of Blackwell s 
f^rpinriorta . that the composite of the sensor itself provides no resistance to the heat. Secondly, die 
anmicability of his “combined” equations that account for the contact resistivity are questioned ill this 

in orfer to gain better information about the contact resistivity, this work rederives Blackwell s 
short time equation and indudes the heat propagation into the sensor. Tnis equation attempts to provide 
more information about the time regime where contact resistivity is d omin a n t in order to eliminate 

the observed variance. 


fimwmMvnf Results 

The preliminary resulting short time equation with first term approximations results in the following 
fqimtinn in tfic I api**** domain: 


T 

P Tj 


i Bio -1/2 


1 • Bio -- 1/2 
1+ t7*’ 



( 1 ) 


1 Blackwell, JJL A Transleni-Flow Method for Determination of Thermal Constants of Insulating Materials in Bulk. JfflirM l flf 
A pplied Ptrwica. Vol. 25, Number 2, Fehnmry, 1954. p. 137-144. 
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and Q’ is the power supplied per unit area of the sensor, R is the radius of die sensor, Xi is die thermal 

conductivity of the sensor, T c is the characteristic temperature, H is the heat tran^er coemdait, Xj is the 
SS conductivity of the medium, D t is the diffusivity of the probe and Dj « > the difTmvity of the 
ryi^Untn Equation (1) represents a more accurate approximation than Bl ackwe ll s equation. If the third 
term of d*e3enominatOTmBquation (1) is dropped, the denominator u 

CT pnnfl"", t m d Anally the inverse Laplace transform is taken the equation reduces to the exact equation 
Blackwell presented: 
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The next step was to take the second term approximation which resulted in the following equation 
which is still in the Laplace transform domain: 
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shortSme temperature profile 

solutionpr6sec^s?by BlackwelL Although the time dependence remains 

aignriicamd^atitm This deviation provides pramse feat a wU1 obtamcd ^ 

ultimately lead to an isolation of the variance provided by die contact resistivity. 


rw«m11 Summary 

formi improved mediums. This inconsistency in data 

acouStion has'^^hypo^^S^to be due to a variance in contact rcastivityat die probe’s surface because 


Incorporating tin 
s work. The fust 
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with different coefficients. The continuing wcrkSSSdp^SSa more accurate mathematical 
^hKJSl^^Se isolation of contact resistivity’s influence and thus allow a consistent and 
H /y»nrnfr> thermal conductivity to be evaluated. 
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MOISTURE MEASUREMENTS 

The most important indication of water state (especially m microgravity) would be the direct 
measurementofthe free energy status of the water in the system. Given the free c .^Syd ls ^uLioii of 
the water in a root module, movement and availability to plants could be easily calculated. Matrix 
nntential ( energy level) can be measured directly by a tensiometer, but these sensors are typically bulky, 

L spice ye. (Ref. d). Since moisture content end 

matrix poSatial are functions of one another, matrix potential is usually inferred from water content 
measurements. 

ril __, ««t*-nf-th£-Hrt methods for determining sub-strate water content include neutron attenuation, 

S^Time Domein Reflectomeuy (TOR). end tarn 
twome^odeTneutron iSJuetion end gemme r.y ittenuetion. pose . mdmtion haza d a nd read to be 
v n jt fv devices Gvosum blocks have been available for ter rest rial measurements for years* bu they 

Wab freflUMicies and are not overly sensitive, The TDR method has become popular because It 
SSS^SesmidW Traditional TOR does require complex circuits, however, 
is nexibie. n . 3 GHz) which would cause concern aboard a spacecraft. Smaller, 

u^u a^ S ^veloped which may Increase their potential for 
J^Tum bfock method relates elecScal to wro - content, the hem pulse 
method uses changes in thermal conductivity to infer moisture level (Ref. 4). 

TOR^tTOCCTStim ^i™ 1 ?mvestipMd! ^i^wMde^S'thaftte con and long qualification time 
rrn uinsTfo/acceptince of this type sensor made it unacceptable. Limitation m power requirements, 

.board a spwtecreft pUcc.smngent Mi « • «**»» 

Mid^reoundingtM^, marred to as contact resistance, can affect a heat pulse measurement of 
moisture level within a subsume. 

n heat nulse sensors have been developed with a solid porous ceramic coating (Ref. 6 ). These 

wau^coment of tte ceramic coating witofo «n»r temperature me -B ecame to* 
SmaiVatilrT^of the sensor is dependent on the water content of the ceramic coaung, other factors like 

impact on the sensor readings. When placed in a soil, the water potential of 
C^KSSettoS ^iSnU the free energy level of toe wrner in. toe soil, H pee. dre ymaor 
rv«m he calibrated directlv as a function of the water potential of the surrounding media. Because the 
SStaiIhlV«utm?Sd rigid contact with die ceramic, the problems associate! w.th toe use of bare 
probes in varying sized media (contact resistance) Is eliminated. 

Airhrtitffh the ceramic surrounding the probe has a known temperature rise to water potential (water 

^d‘^" The°prc^e ^iththc 
caS^tc TW is Dossible only if the hydraulic conductivity of the small grained ceramic has a simitar 
^SSS^SiSS^SSS^mZ substrate. The substrate u«d *» 

requirement. Selecting smaller substrates will eventually allow the direct measurement of matrix 

potential using this sensor. 
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Propagation of menisci into pore space: 
Simulation by finite elements 

Jean Ro 


The finite element method has been successfully applied to non-linear fluid mechanics problems. 
Finite element algorithms make use of the weak forms of the basic equations of motion and 
continuity equations. The solution is obtained as an approximation to the unknown functions 
using trial functions defined over individual elements. The main advantage of the finite element 
method is that the elements do not have to be of uniform size or shape. Therefore, irregular 
domains require no special handling - a particular advantage in dealing with nee boundary 
problems such as ours. In addition, more efficient computation can be achieved by using smaller 
elements in regions of rapid variations of the dependent variables, or by raising the degree of the 
polynomials which are used as trial functions. 

For application to free boundary problems, or when the fluid is viscoelastic, it is useful to work 
with stress as a dependent variable. The weak form of the equation of motion contains stresses 
which are implicit functions of the velocity field, since in the general non-Newtonian case, stress 
components cannot be given as algebraic functions of the velocity components and their 
derivatives. Prof. Khomami's finite element code treats the deviatonc stress tensor as the sum of 
the polymeric (Jp) and Newtonian (T*) solvent stress tensor with the two related to 
appropriate flow variables by a given constitutive relations. In the mixed method, the extra 
stresses are treated as unknowns, and are assumed as linear combinations of tnal functions just as 
the velocity and the pressure are treated. Newtonian flow problems can be treated simply by 
setting the polymeric stress components to zero and using the Newtoman relatio^iii tetw^n 
stress and rate of strain. For viscoelastic flow problems, this mixed method is combined with the 
elastic-viscous split stress technique (EVSS) in which the deviatonc stress is divid^ntoelasuc 
and viscous portions, with the viscous part including the Newtonian conmbuuoi^ tom twth the 
solvent part of the stress and the polymer part. The advantage of adopting the EVSS method is that 
convergent solutions can be obtained regardless of the value of the ratio between the so *^ nl . 
viscosity and the viscosity of the solution since equations are rewritten in an explicitly elliptic way. 
The overall strategy in the code is to solve the continuity and momentum equations, which contain 
an elliptic operator, by Galeridn's method, and to compute stresses by Streamline Unwinding 
Petrov-Galeririn method (SUPG) which is suitable for hyperbolic equations. 
elements arc allowed. Biquadratic shape functions are used for stresses and velocity, and bilinear 
shape functions are used for the rate of strains and pressure. There are total of ten variables per 

element. 

The finite element method is well suited for our free surface problem with the irregular domain. In 
our iterative formulation, the finite element solver mentioned above is used to soiveforthe two 
dimensional flow with the kinematic and tangential stress boundary conditions imposed on the 
ftgc mnftH free boundary which is unknown a-priori. Then the computed pressure and normal 
stresses together with the normal stress balance, are used to corrcctthe free boundary location. 
That is accomplished by solving the differential equation resulting from the normal swss balance 
equation using the finite element method with one dimensional free surface as our domain. The 
normal stress condition depends on the curvature of the free surface, whereas |*ebn^«cand 
tangential stress boundary conditions depend only on surface orientation. Therefore, the nonnal 
stress scheme mentioned above is more likely to converge than other possible schemes. 

We are currently in the process of validation of this procedure by solving the problem of meniscus 
propagation between parallel plates at finite capillary number. This will be followed by a study of 
meniscus propagation between spatially varying walls more appropriate to geometric models of 

porous media. 



Flow of Water through a Bed of Packed Spheres 

A. Farooq 

Dept, of Mechanical Engineering, Stanford University, Stanford, CA 94305. 


1 Introduction 

Some recent experiments (Yendler et al. 1993) have investigated the propagation of water 
through a bed of packed spheres. The motivation is to investigate the propagation of water in 
micro- gravity environments with the objective of growing plants in space. Their experimental 
set up is shown schematically in Figure 1. It consists of a long channel of length L and height 
H which is packed with spheres of diameter D p . This packed bed is fed with water (through 
a mechanism not shown, see original paper for more details) without any hydraulic head, so 
that the driving force for propagation is dominated by capillarity. The experiment consists 
of recording the movement of water as it propagates through an initially dry bed. 

The geometric parameters of the experiment are the dimensions of the set up, (given as 
L = 188.9mm, H = 12.7mm, from Yendler, 1993), the sphere diameter, D P for which values 
of 0.5, 0.75, 1.0, 1.5, 2.0 (all diameters in mm), and the porosity of the bed e which has been 
reported to be in the range e « 0.395 - 0.421. Note that in any given experiment, beads of 
only a single size are used. 

Defining the saturation 9, to be the fractional water content in the bed, ($mo* = e )> the 
experiment records S which is some suitable value of 9 , (say 0.99e) as a function of time 
t. As is well known from the theory of two-phase immiscible capillary flow through porous 
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media (to be developed in the next section), 8 ~ \/i. 

In figure 2 we have redisplayed the results obtained by Yendler (1993). Figure 2a shows a 
plot of 8 as a function of t for smaller sphere sizes, (D p = 0.5, 0.7, 1.0) and 2b shows the 
same for D p = 1.5, 2.0, mm. The interesting observation is that the rate of propagation is 
lower for the larger spheres. In the next section we will develop the equation governing the 
flow of liquid in unsaturated media and try to understand this peculiar behaviour. 


2 Mathematical Formulation of Governing Equations 

The equation governing the propagation of water through packed beds is the Richard’s 
equation, given below (in dimensional form); 


86 

dt 


V.DV9 + 


(pg\ dkdl 

\fi)d6dy 


( 1 ) 


where D = k$. Here if) is the capillary pressure and is a function of 6 alone, k is the 
permeability of the medium. Hence the diffusion coefficient D, is also a function of 9 making 
(1) a non-linear diffusion equation. The second term on the right hand side of (1) accounts 
for the effect of buoyancy. Richard’s equation is based on the principle of continuity and 
D’Arcy’s law. 1 

iThe continuity equation is ff + V.u = 0 and D’Arcy’s law is given by u = - V<f>, for a suitable 

pressure <j>, which in our case can be decomposed into capillarity and hydraulic heads: <(> = tp(0) + pgy. 
Combining the above three expressions yields (1) after some algebraic rearrangement. 
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According to Buckley-Leverett theory (see Bear 1972), the capillary head can be related to 
the surface tension and permeability of the medium as follows: 


* - ( 2 ) 

where <t is the surface tension and J{6) is a function characteristic of the medium. The 
permability k , is a function of both the nature of the porous media as well as the saturation 
6. It is common to scale the permeability with respect to the characteristic grain size of the 
medium, (in our case, D p ) and a relative permeability to account for differences in saturation 
as, 

Vk = D P G{6 ) (3) 

where G(9) is the relative permeability. We note that both G(6) and J(6) are 0(1). 

2.1 Non-dimensionalization 

Choosing the characteristic length scale as H , the height of the apparatus, and as the 
characteristic time scale, (1) cam be cast as: 

= V.GJ'VB + 2 GG'Bo^ (4) 

where r is the rescaled time, amd X, Y are the rescaled coordinates. The dimensionless 
group, Bo = paB a D r is the Bond number and gives the relative importance of gravity and 
surface tension (capillarity). We note that D has the scale and thus the diffusivity 

scales linearly with D p . 
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2.2 Flow with Bo — * 0 


Setting V = G(9)J'(9) in (4) we obtain in the limit Bo — ► 0 


-S- = V.PV0 
or 


Seeking a one dimensional solution to (5), we write it as: 


86 _ 9 v 89 

dr ~ dx v dx 


with inital and boundary conditions given by 


( 5 ) 


( 6 ) 


9 = 0 for X > 0, t = 0 
9 = e for X = 0, r > 0 


Equation (5) with associated initial and boundary conditions admits a similarity solution in 
the similarity variable, rj = ^ and 


d ^df 7i df n 

— V— + = 0 

drf dr} 2 drj 


(7) 


where f(r]) = 9 (X,t). Solutions to this ODE can be obtained after making suitable assump- 
tions for the functional dependence of V = V{9). If it is assumed that V is a constant, then 
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D p) mm 

m 

IS 

6 

11 

D 

0.46 

1.194 

1.425 

3.097 

0.75 

1.533 

2.350 

3.130 

1.00 

1.493 

2.229 

2.229 

■1 

1.50 

0.723 

0.273 

0.182 

2.00 

0.100 

0.010 

0.005 


Table 1: This table shows the effective diffusivity, D and examines the validity of its scaling 
with respect to D p . 

the solution has the form of the familiar complimentary error function. For V = D(6), the 
equation becomes non-linear and generally needs to be solved numerically. This has been 
done by P hilli p (1955), which can be consulted for more details. The results in Phillip (1955) 
show that the basic structure of the solutions remains the same as for linear diffusion, the 
exact profile being a strong function of the diffusion coefficient. Hence the locus of points for 
£ 0 = 0.99 behaves as as ~ as has been observed in Figure 2(a,b), where the propagation 

rate is shown to scale as y/t 2 . Yendler et al (1993) also calculate the effective diffusivity, 
V by fitting a line l = rny/t + b to the data of figure 2. Their results are given in table (1) 
where D = m 2 , is the effective (dimensional) diffusivity. 

The last col umn of the table examines the scaling of the effective diffusivity D with respect 
to D p , the grain size and surprisingly, the relationship does not hold for D p > 0.75. The 
only assumption in equation (1) is is D’Arcy’s law which assumes that the Reynolds number 

2 Yendler et al have presented their results for unsealed (dimensional time) 
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D p , mm 

U (cm/min) 

Ren, 

0.50 

33.33 

2.75 

0.75 

20.00 

2.50 

1.50 

6.00 

1.50 

2.00 

0.26 

0.53 


Table 2: This table shows the Reynold’s number, ReD r for the range of D p considered. 

of the flow is low. There is considerable experimental data (see for example Bird, Stewart, 
Lightfoot , 1960, page 198) which suggests that D’Arcy’s law is valid for a pore Reynolds 
number, R^d, < 10. 

Using data available from the experiments, we have estimated Rev r — using v = 

10 -6 m 2 s -1 and found the value to be less than 3 for all the pore sizes considered from 0.50 
to 2.0. This is shown in table (2) where the velocity (as obtained from experiments) and the 
calculated Reynolds number are given 

The capillary number, Ca = has been found to be of the order of 10 -4 (using cr = 

70 X 10 -3 Nm -1 , and a reference velocity scale). Thus the capillary forces are dominant 
with respect to at the pore scales. Hence we conclude that (1) should correctly predict the 
behaviour of the fluid in the given range of sphere diameters. 
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D p , mm 

Bo 

0.50 

1.57 

0.75 

2.35 

1.00 

3.14 

1.50 

4.71 

2.00 

6.28 


Table 3: This table shows the Bond number, Bo for the range of D p considered. 

2.3 Effect of Bond number, Bo 


Horizontal propagation of water through a bed of packed spheres was studied with the 
implicit assumption that gravity is not an important factor, and hence the model could be 
used to study similar situations in microgTavity environments. This assumption now needs 
to be re-evaluated in light of the discrepency observed between one- dimensional theory as 
developed in the previous section and the experimental results. 

Using p = 10 3 , g = 10, a = 70 x 10" 3 , and H = 12 x 1CT 3 , (all in SI units) the Bond number 
is calculated and shown in table (3). We note that Bo > 2.5 for the larger spheres and hence 
gravity is significant in (1) and cannot be ignored. 
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3 Conclusion and suggestion for further work 

We have shown that the Bond number is significantly large (greater than 1) for the larger 
sphere diameters under consideration and hence gravity plays a significant role in the outcome 
of the experiment. This could be responsible for the diminished rate of propagation that has 
been observed. 3 

It may be possible to overcome this difficulty, simply by controlling the Bond number. Since 
we may desire to maintain the grain size and the size of the apparatus at current values, the 
only parameters left are the liquid density and surface tension. Hence it may be possible 
to conduct the experiment successfully by a suitable choice of liquid, with higher surface 
tension and lower density. Table 3 suggests that a factor of 6 — 10 reduction in Bond number 
would be sufficient. 
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Figure 2: Experimental results obtained by Yendler (1993). 
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